Se tiene dataset consitutido por observaciones de 11 variables en 34 estados del mundo. Las variables se encuentran estandarizadas pues están tomadas con unidades de medida muy diferentes. Las variables son:
Mostramos los datos:
library(foreign)
data <- read.spss("./DB_3.sav", to.data.frame = TRUE)
data
Comprobamos que las variables están normalizadas:
summary(data)
## PAIS ZPOBDENS ZTMINFAN ZESPVIDA
## Length:34 Min. :-1.0778 Min. :-1.1026 Min. :-2.1453
## Class :character 1st Qu.:-0.8497 1st Qu.:-0.9586 1st Qu.:-0.7460
## Mode :character Median :-0.1616 Median :-0.3931 Median : 0.2781
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5547 3rd Qu.: 0.8985 3rd Qu.: 0.9033
## Max. : 2.8616 Max. : 1.9048 Max. : 1.2486
##
## ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## Min. :-1.7697 Min. :-1.1473 Min. :-1.2342 Min. :-1.88521
## 1st Qu.:-0.7320 1st Qu.:-0.8829 1st Qu.:-0.8480 1st Qu.:-0.72858
## Median : 0.1268 Median :-0.2916 Median :-0.2134 Median : 0.03541
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.8148 3rd Qu.: 0.8694 3rd Qu.: 0.7115 3rd Qu.: 0.85176
## Max. : 1.5096 Max. : 2.3717 Max. : 1.9052 Max. : 1.62885
##
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## Min. :-0.9696 Min. :-0.86586 Min. :-2.1341 Min. :-0.9507
## 1st Qu.:-0.9240 1st Qu.:-0.59889 1st Qu.:-0.6735 1st Qu.:-0.7813
## Median :-0.3237 Median :-0.20626 Median :-0.1067 Median :-0.3900
## Mean : 0.0000 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.7736 3rd Qu.: 0.07996 3rd Qu.: 0.8789 3rd Qu.: 0.5828
## Max. : 2.4024 Max. : 4.42620 Max. : 1.7045 Max. : 2.7498
## NA's :1
Mostramos desviación típica de las variables
print("Desviación típica:")
## [1] "Desviación típica:"
sapply(data, sd)
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm =
## na.rm): NAs introduced by coercion
## PAIS ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## NA 1 1 1 1 1 1 1
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## NA 1 1 1
Desviación típica de la variable ZTLIBROP sin el valor perdido:
print("Desviación típica:")
## [1] "Desviación típica:"
sd(data[,9][-22])
## [1] 1
Recodificamos a UTF-8 y aprovechamos para corregir el nombre de algunos paises incompletos y/o con erratas.
data[1, "PAIS"] <- "áfrica sur"
data[3, "PAIS"] <- "argentina"
data[4, "PAIS"] <- "australia"
data[11, "PAIS"]<- "españa"
data[14, "PAIS"]<- "hungría"
data[16, "PAIS"]<- "indonesia"
data[20, "PAIS"]<- "japón"
data[22, "PAIS"]<- "marruecos"
data[23, "PAIS"] <- "méxico"
data[27, "PAIS"] <- "rd alemania"
data[28, "PAIS"] <- "reino unido"
data[29, "PAIS"] <- "rf alemania"
Vemos los valores perdidos:
cbind(apply(is.na(data),2,sum),apply(is.na(data),2,sum)/dim(data)[1])
## [,1] [,2]
## PAIS 0 0.00000000
## ZPOBDENS 0 0.00000000
## ZTMINFAN 0 0.00000000
## ZESPVIDA 0 0.00000000
## ZPOBURB 0 0.00000000
## ZTMEDICO 0 0.00000000
## ZPAGRICU 0 0.00000000
## ZPSERVI 0 0.00000000
## ZTLIBROP 1 0.02941176
## ZTEJERCI 0 0.00000000
## ZTPOBACT 0 0.00000000
## ZTENERGI 0 0.00000000
vemos un valor perdido en ZTLIBROP. Aprovechamos para guardar una copia del dataset sin variable PAIS que no es de interés.
data_ <- data[-1]
head(data_)
Los datos perdidos en ZTLIBROP son inferiores al 5%, sustituiremos por la media:
not_available<-function(data,na.rm=F){
data[is.na(data)]<-mean(data,na.rm=T)
data
}
data_<-as.data.frame(apply(data_, 2, not_available))
Veamos un resumen de las variables:
summary(data_)
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB
## Min. :-1.0778 Min. :-1.1026 Min. :-2.1453 Min. :-1.7697
## 1st Qu.:-0.8497 1st Qu.:-0.9586 1st Qu.:-0.7460 1st Qu.:-0.7320
## Median :-0.1616 Median :-0.3931 Median : 0.2781 Median : 0.1268
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5547 3rd Qu.: 0.8985 3rd Qu.: 0.9033 3rd Qu.: 0.8148
## Max. : 2.8616 Max. : 1.9048 Max. : 1.2486 Max. : 1.5096
## ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## Min. :-1.1473 Min. :-1.2342 Min. :-1.88521 Min. :-0.9696
## 1st Qu.:-0.8829 1st Qu.:-0.8480 1st Qu.:-0.72858 1st Qu.:-0.9145
## Median :-0.2916 Median :-0.2134 Median : 0.03541 Median :-0.2442
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.8694 3rd Qu.: 0.7115 3rd Qu.: 0.85176 3rd Qu.: 0.6923
## Max. : 2.3717 Max. : 1.9052 Max. : 1.62885 Max. : 2.4024
## ZTEJERCI ZTPOBACT ZTENERGI
## Min. :-0.86586 Min. :-2.1341 Min. :-0.9507
## 1st Qu.:-0.59889 1st Qu.:-0.6735 1st Qu.:-0.7813
## Median :-0.20626 Median :-0.1067 Median :-0.3900
## Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.07996 3rd Qu.: 0.8789 3rd Qu.: 0.5828
## Max. : 4.42620 Max. : 1.7045 Max. : 2.7498
Mostramos boxplot de las variables para visualizar su distribución:
boxplot (data_, main = "Boxplot de las variables", xlab = "Variables", ylab = "Valores")
Mostramos coef. de asimetría de las variables
library(moments)
skewness(data_)
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## 1.0679539 0.5477565 -0.5802673 -0.3189850 0.4909934 0.5918253 -0.1239223
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## 0.8170779 2.8653904 -0.1284517 1.2270338
Mostramos curtosis de las variables
kurtosis(data_)
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## 3.542159 1.796001 2.156273 1.902481 2.065480 2.055664 1.911288 2.524152
## ZTEJERCI ZTPOBACT ZTENERGI
## 12.634918 2.111551 3.912472
Analizaremos con más detalle:
ZPOBDENS
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 1.06373e-16
##
##
## $centralidad$resistente
## [,1]
## PMC -0.1474835
## Trimedia -0.1545405
## Centrimedia -0.2293829
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 9.400883e+15
##
## $dispersion$clasica$rango
## [1] -1.077834 2.861621
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.4043955
## MEDA 0.6924864
## CVc -4.7611940
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 1.067954
##
## $forma$clasica$kurtosis
## [1] 3.542159
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.08733974
## Asimetría de Kelly -0.37369403
## Asimetría de Kelly adimensional 2.31250000
Las medidas resistentes de centralidad están ligeramente hacia la izquierda. MEDA inferior a desviación típica. Observamos asimetría, así como cierta acumulación de datos por el valor de curtosis.
ZTMINFAN
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 8.982296e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.0300511
## Trimedia -0.2115862
## Centrimedia -0.2268480
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 1.113301e+16
##
## $dispersion$clasica$rango
## [1] -1.102554 1.904824
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.8571199
## MEDA 0.6040459
## CVc -30.8993711
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.5477565
##
## $forma$clasica$kurtosis
## [1] 1.796001
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.9235577
## Asimetría de Kelly -0.6132994
## Asimetría de Kelly adimensional 1.5600769
Observamos cierto desplazamiento a la izquierda en las medidas resistentes de centralidad. Asimetría algo menor que el caso anterior, así como una distribución más aplanada.
ZESPVIDA
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -1.72884e-16
##
##
## $centralidad$resistente
## [,1]
## PMC 0.07863105
## Trimedia 0.17836465
## Centrimedia 0.17613181
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] -5.784225e+15
##
## $dispersion$clasica$rango
## [1] -2.145279 1.248640
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.6493257
## MEDA 0.7621433
## CVc 10.4877506
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.5802673
##
## $forma$clasica$kurtosis
## [1] 2.156273
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.7172544
## Asimetría de Kelly 0.4995611
## Asimetría de Kelly adimensional 1.7963476
Observamos distrbución aplanada aunque menos que la anterior, variación en los datos, así como medidas resistentes de centralidad ligeramente desplazadas a la derecha.
ZPOBURB
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 4.511206e-16
##
##
## $centralidad$resistente
## [,1]
## PMC 0.0413792
## Trimedia 0.0840905
## Centrimedia 0.1147624
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 2.216702e+15
##
## $dispersion$clasica$rango
## [1] -1.769694 1.509616
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5467796
## MEDA 0.7223655
## CVc 18.6903015
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.318985
##
## $forma$clasica$kurtosis
## [1] 1.902481
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.6736702
## Asimetría de Kelly 0.2958259
## Asimetría de Kelly adimensional 2.3329787
Observamos distribución con cierta asimetría, la más aplanada tras la de ZTMINFAN, y prácticamente no hay desplazamiento de centralidad.
ZTMEDICO
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 2.416495e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.006716126
## Trimedia -0.149133349
## Centrimedia -0.149734266
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 4.138224e+16
##
## $dispersion$clasica$rango
## [1] -1.147256 2.371712
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.7522727
## MEDA 0.7980172
## CVc -130.4526316
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.4909934
##
## $forma$clasica$kurtosis
## [1] 2.06548
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.9769641
## Asimetría de Kelly -0.3735297
## Asimetría de Kelly adimensional 1.2811833
Se aprecia distribución aplanada, cierta asimetría, así como ligero desplazamiento a la izquierda.
ZPAGRICU
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 6.355338e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.06825485
## Trimedia -0.14081091
## Centrimedia -0.20433276
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 1.57348e+16
##
## $dispersion$clasica$rango
## [1] -1.234234 1.905157
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5595319
## MEDA 0.7069276
## CVc -11.4243309
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.5918253
##
## $forma$clasica$kurtosis
## [1] 2.055664
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.6801059
## Asimetría de Kelly -0.4817497
## Asimetría de Kelly adimensional 2.2578456
Distribución aplanada, cierta asimetría, así como medidas resistentes de centralidad ligeramente desplazadas a la izquierda.
ZPSERVI
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -3.529078e-17
##
##
## $centralidad$resistente
## [,1]
## PMC 0.0615893232
## Trimedia 0.0485015920
## Centrimedia 0.0006495749
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] -2.833601e+16
##
## $dispersion$clasica$rango
## [1] -1.885211 1.628845
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5803435
## MEDA 0.7918077
## CVc 12.8296875
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.1239223
##
## $forma$clasica$kurtosis
## [1] 1.911288
##
##
## $forma$resistente
## 25%
## Asimetría de Yule 0.73913043
## Asimetría de Kelly 0.05071496
## Asimetría de Kelly adimensional 1.43206522
Distribución aplanada, prácticamente centrada, y poca asimetría.
ZTLIBROP
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 6.729069e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.1111009
## Trimedia -0.1776580
## Centrimedia -0.2615358
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 0.9847319
##
## $dispersion$clasica$Coef_varización
## [1] 1.4634e+16
##
## $dispersion$clasica$rango
## [1] -0.9696156 2.4023604
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.6067690
## MEDA 0.6849277
## CVc -7.2311230
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.8170779
##
## $forma$clasica$kurtosis
## [1] 2.524152
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.5450695
## Asimetría de Kelly -0.5339910
## Asimetría de Kelly adimensional 2.1865595
Se observa fuerte asimetría hacia la derecha, así como desplazamiento a la derecha y concentración de los datos por el alto valor de curtosis. ZTEJERCI
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 1.859184e-16
##
##
## $centralidad$resistente
## [,1]
## PMC -0.2594621
## Trimedia -0.2328606
## Centrimedia -0.2192510
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 5.378705e+15
##
## $dispersion$clasica$rango
## [1] -0.8658606 4.4262018
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 0.6788462
## MEDA 0.3313997
## CVc -1.3081800
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 2.86539
##
## $forma$clasica$kurtosis
## [1] 12.63492
##
##
## $forma$resistente
## 25%
## Asimetría de Yule 0.2579423
## Asimetría de Kelly -0.3448882
## Asimetría de Kelly adimensional 1.6721112
Observamos muy fuerte asímetría a la derecha, desplazamiento de las medidas resistentes de centralidad a la izquierda, así como alta concentración de los datos por valor muy alto de curtosis. ZTPOBACT
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -2.952027e-16
##
##
## $centralidad$resistente
## [,1]
## PMC 0.10268877
## Trimedia -0.00198850
## Centrimedia 0.02917055
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] -3.387502e+15
##
## $dispersion$clasica$rango
## [1] -2.134094 1.704472
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5523378
## MEDA 0.8557514
## CVc 7.5584589
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.1284517
##
## $forma$clasica$kurtosis
## [1] 2.111551
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -1.96271531
## Asimetría de Kelly -0.06440751
## Asimetría de Kelly adimensional 0.60382547
Distribución aplanada, poca asimetría y prácticamente centrada. ZTENERGI
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 3.586955e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.09924855
## Trimedia -0.24461072
## Centrimedia -0.26271176
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 2.78788e+16
##
## $dispersion$clasica$rango
## [1] -0.950661 2.749785
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.364062
## MEDA 0.520457
## CVc -6.871948
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 1.227034
##
## $forma$clasica$kurtosis
## [1] 3.912472
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.7454988
## Asimetría de Kelly -0.5056537
## Asimetría de Kelly adimensional 1.2966382
Se observa distribución algo desplazada a la izquierda, concentrada (alto valor de curtosis), y asimétrica a la derecha.
Veamos de nuevo boxplot de las variables:
boxplot (data_, main = "Boxplot de las variables", xlab = "Variables", ylab = "Valores")
En cierta manera todas las variables parecen tener cierta similitud a una normal. Vemos 1 outlier en ZPOBDENS y ZTENERGI, varios en ZTEJERCI. Los sustituiremos por la media.
outlier<-function(data,na.rm=T){
H<-1.5*IQR(data)
if(any(data<=(quantile(data,0.25,na.rm = T)-H))){
data[data<=quantile(data,0.25,na.rm = T)-H]<-NA
data[is.na(data)]<-mean(data,na.rm=T)
data<-outlier(data)}
if(any(data>=(quantile(data,0.75, na.rm = T)+H))){
data[data>=quantile(data,0.75, na.rm = T)+H]<-NA
data[is.na(data)]<-mean(data,na.rm=T)
data<-outlier(data)
}
return(data)
}
data_<-apply(data_, 2, outlier)
boxplot(data_, main = "Boxplot de las variables", xlab = "Variables", ylab = "Valores")
Mostramos resumen de los datos tras tratamiento de outliers:
summary(data_)
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB
## Min. :-1.07783 Min. :-1.1026 Min. :-2.1453 Min. :-1.7697
## 1st Qu.:-0.84968 1st Qu.:-0.9586 1st Qu.:-0.7460 1st Qu.:-0.7320
## Median :-0.16160 Median :-0.3931 Median : 0.2781 Median : 0.1268
## Mean :-0.08672 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.40244 3rd Qu.: 0.8985 3rd Qu.: 0.9033 3rd Qu.: 0.8148
## Max. : 2.15204 Max. : 1.9048 Max. : 1.2486 Max. : 1.5096
## ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## Min. :-1.1473 Min. :-1.2342 Min. :-1.88521 Min. :-0.9696
## 1st Qu.:-0.8829 1st Qu.:-0.8480 1st Qu.:-0.72858 1st Qu.:-0.9145
## Median :-0.2916 Median :-0.2134 Median : 0.03541 Median :-0.2442
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.8694 3rd Qu.: 0.7115 3rd Qu.: 0.85176 3rd Qu.: 0.6923
## Max. : 2.3717 Max. : 1.9052 Max. : 1.62885 Max. : 2.4024
## ZTEJERCI ZTPOBACT ZTENERGI
## Min. :-0.86586 Min. :-2.1341 Min. :-0.9507
## 1st Qu.:-0.59889 1st Qu.:-0.6735 1st Qu.:-0.7813
## Median :-0.28969 Median :-0.1067 Median :-0.3900
## Mean :-0.28969 Mean : 0.0000 Mean :-0.1690
## 3rd Qu.:-0.04621 3rd Qu.: 0.8789 3rd Qu.: 0.3346
## Max. : 0.68708 Max. : 1.7045 Max. : 1.5672
y los datos iniciales:
summary(data)
## PAIS ZPOBDENS ZTMINFAN ZESPVIDA
## Length:34 Min. :-1.0778 Min. :-1.1026 Min. :-2.1453
## Class :character 1st Qu.:-0.8497 1st Qu.:-0.9586 1st Qu.:-0.7460
## Mode :character Median :-0.1616 Median :-0.3931 Median : 0.2781
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5547 3rd Qu.: 0.8985 3rd Qu.: 0.9033
## Max. : 2.8616 Max. : 1.9048 Max. : 1.2486
##
## ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## Min. :-1.7697 Min. :-1.1473 Min. :-1.2342 Min. :-1.88521
## 1st Qu.:-0.7320 1st Qu.:-0.8829 1st Qu.:-0.8480 1st Qu.:-0.72858
## Median : 0.1268 Median :-0.2916 Median :-0.2134 Median : 0.03541
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.8148 3rd Qu.: 0.8694 3rd Qu.: 0.7115 3rd Qu.: 0.85176
## Max. : 1.5096 Max. : 2.3717 Max. : 1.9052 Max. : 1.62885
##
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## Min. :-0.9696 Min. :-0.86586 Min. :-2.1341 Min. :-0.9507
## 1st Qu.:-0.9240 1st Qu.:-0.59889 1st Qu.:-0.6735 1st Qu.:-0.7813
## Median :-0.3237 Median :-0.20626 Median :-0.1067 Median :-0.3900
## Mean : 0.0000 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.7736 3rd Qu.: 0.07996 3rd Qu.: 0.8789 3rd Qu.: 0.5828
## Max. : 2.4024 Max. : 4.42620 Max. : 1.7045 Max. : 2.7498
## NA's :1
Veamos si estamos tratando con variables normales, para ello utilizamos el método gráfico qq-plot.
par(mar=c(1,1,1,1))
par(mfrow=c(3,5))
invisible(apply(data_, 2, function(x){
qqnorm(x,main=NULL)
abline(a=0,b=1,col="red")
}))
Vemos que las variables que más se acercan a la normalidad son la 4, 6, 7 y 10 (ZPOBURB, ZPAGRICU, ZPSERVI y ZTPOBACT), en menor medida la 1, 2, 3, 5, 8 y 11 (ZPOBDENS, ZTMINFAN, ZESPVIDA, ZTMEDICO, ZTLIBROP y ZTENERGI), la que menos se acerca a normalidad es la 9 (ZTEJERCI).
Normalizamos los datos
data_norm <- scale(data_)
Trabajaremos con los datos normalizados:
datos_pca <- data_norm
row.names(datos_pca) = data$PAIS
Mostramos las correlaciones entre variables
library(corrr)
datos_pca_correlaciones <- correlate(datos_pca) #Cálculo de objeto de correlaciones
##
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'
rplot(datos_pca_correlaciones, legend = TRUE, colours = c("firebrick1", "black","darkcyan"), print_cor = TRUE)
## Don't know how to automatically pick scale for object of type noquote. Defaulting to continuous.
## Don't know how to automatically pick scale for object of type noquote. Defaulting to continuous.
##
## Attaching package: 'psych'
## The following object is masked _by_ '.GlobalEnv':
##
## outlier
Hacemos test de esfericidad de Bartlett por el que comprobaremos si las correlaciones son significativamente distintas de 0
cortest.bartlett(cor(datos_pca), nrow(datos_pca))
## $chisq
## [1] 387.5628
##
## $p.value
## [1] 1.583951e-51
##
## $df
## [1] 55
Observamos un p-valor muy bajo, se rechaza hipótesis nula por lo que los datos no están incorrelados.
Haremos Análisis de Componentes Principales, los datos ya los hemos preparado para ello.
PCA<-prcomp(datos_pca)
Mostramos los coeficientes de las componentes principales, esto es, el peso de cada variable en la correspondiente componente principal:
PCA$rotation
## PC1 PC2 PC3 PC4 PC5
## ZPOBDENS 0.06796131 0.39036079 0.07578571 -0.873056137 0.05981546
## ZTMINFAN -0.37529798 -0.11573027 -0.08212636 -0.014741875 0.41788865
## ZESPVIDA 0.37165513 0.07270402 -0.00407411 0.012970861 -0.53281544
## ZPOBURB 0.35808785 -0.29124050 -0.06094609 -0.113795410 -0.05152603
## ZTMEDICO 0.33374417 0.13909259 -0.09587119 0.280930962 0.15782496
## ZPAGRICU -0.36224935 0.29940891 0.04348437 0.061235615 -0.15998056
## ZPSERVI 0.29607632 -0.48878963 0.09183567 -0.148420697 -0.03064997
## ZTLIBROP 0.32594252 0.08892433 -0.07078859 -0.135187927 0.28062939
## ZTEJERCI 0.02397325 0.17948476 -0.93319975 -0.006241743 -0.12518016
## ZTPOBACT 0.20289865 0.56167354 0.29659750 0.281932179 -0.11796713
## ZTENERGI 0.33155119 0.20152127 -0.02054058 0.148110844 0.61275026
## PC6 PC7 PC8 PC9 PC10
## ZPOBDENS -0.24337351 0.09619809 -0.04353317 0.022846323 -0.01646372
## ZTMINFAN -0.01756878 0.05245475 -0.39097670 -0.157260993 -0.47091956
## ZESPVIDA 0.03028763 -0.04105430 0.29289664 -0.091799724 -0.54354628
## ZPOBURB -0.07065324 -0.10651593 -0.29868358 -0.792450091 0.16228759
## ZTMEDICO -0.36060911 0.78403566 -0.06560314 0.001929436 -0.06125296
## ZPAGRICU 0.13248858 0.09343760 -0.01200167 -0.325592955 -0.48343682
## ZPSERVI -0.07352269 -0.09280074 -0.37422057 0.449245272 -0.39418278
## ZTLIBROP 0.85707761 0.21913126 -0.01229527 0.015266047 -0.02755826
## ZTEJERCI -0.02870986 -0.16141551 -0.18218984 0.127931321 0.01085377
## ZTPOBACT 0.02948798 -0.25103160 -0.61023870 0.097328148 0.08409493
## ZTENERGI -0.21278600 -0.45367225 0.34643517 -0.066593066 -0.23535405
## PC11
## ZPOBDENS -0.016686681
## ZTMINFAN -0.511610948
## ZESPVIDA -0.424556791
## ZPOBURB 0.088144761
## ZTMEDICO 0.058614024
## ZPAGRICU 0.617855535
## ZPSERVI 0.361609762
## ZTLIBROP 0.008370356
## ZTEJERCI 0.063899183
## ZTPOBACT -0.095710436
## ZTENERGI 0.143880098
Mostramos las desviaciones típicas de cada componente principal, proporción de varianza explicada y acumulada:
PCA$sdev
## [1] 2.4818030 1.2805818 1.0311152 0.9515420 0.6448055 0.6027181 0.4818699
## [8] 0.3354192 0.2491626 0.1641974 0.1390789
summary(PCA)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.4818 1.2806 1.03112 0.95154 0.6448 0.60272 0.48187
## Proportion of Variance 0.5599 0.1491 0.09665 0.08231 0.0378 0.03302 0.02111
## Cumulative Proportion 0.5599 0.7090 0.80568 0.88799 0.9258 0.95881 0.97992
## PC8 PC9 PC10 PC11
## Standard deviation 0.33542 0.24916 0.16420 0.13908
## Proportion of Variance 0.01023 0.00564 0.00245 0.00176
## Cumulative Proportion 0.99015 0.99579 0.99824 1.00000
Mostramos gráficamente la varianza explicada
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
varianza_explicada <- PCA$sdev^2 / sum(PCA$sdev^2)
ggplot(data = data.frame(varianza_explicada, pc = 1:11),
aes(x = pc, y = varianza_explicada, fill=varianza_explicada )) +
geom_col(width = 0.3) +
scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
labs(x = "Componente principal", y= " Proporción de varianza explicada")
Y la varianza acumulada
varianza_acum<-cumsum(varianza_explicada)
ggplot( data = data.frame(varianza_acum, pc = 1:11),
aes(x = pc, y = varianza_acum ,fill=varianza_acum )) +
geom_col(width = 0.5) +
scale_y_continuous(limits = c(0,1)) +
theme_bw() +
labs(x = "Componentes principales",
y = "Proporción varianza acumulada")
Para la selección de componentes principales utilizaremos la regla de Abdi et al. (2010). Promediamos las varianzas explicadas por la componentes principales:
PCA$sdev^2
## [1] 6.15934602 1.63988972 1.06319852 0.90543213 0.41577420 0.36326910
## [7] 0.23219856 0.11250601 0.06208202 0.02696079 0.01934294
Calculamos la media:
mean(PCA$sdev^2)
## [1] 1
Y no quedamos con las tres primeras componentes principales pues tienen varianza explicada superior a 1, que es la media.
Visualizaremos ahora la contribución de las componentes principales.
Veamos una comparativa entre la primera y segunda componente principal analizando que variables tienen más peso para la definición de cada componente principal
library("factoextra")
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_pca_var(PCA,
repel=TRUE,col.var="cos2",
legend.title="Distancia")+theme_bw()
Ahora comparativa entre la primera y tercera componente principal analizando que variables tienen más peso para la definición de cada componente principal
fviz_pca_var(PCA,axes=c(1,3),
repel=TRUE,col.var="cos2",
legend.title="Distancia")+theme_bw()
Comparativa entre la segunda y tercera componente principal analizando que variables tienen más peso para la definición de cada componente principal
fviz_pca_var(PCA,axes=c(2,3),
repel=TRUE,col.var="cos2",
legend.title="Distancia")+theme_bw()
Representaremos las observaciones de los objetos junto con las componentes principales, identificando con colores las observaciones que mayor varianza explican de las componentes principales.
Observaciones en la primera y segunda componente principal
fviz_pca_ind(PCA,col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var")+theme_bw()
Observaciones en la primera y tercera componente principal
fviz_pca_ind(PCA,axes=c(1,3),col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var")+theme_bw()
Observaciones en la segunda y tercera componente principal
fviz_pca_ind(PCA,axes=c(2,3),col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var")+theme_bw()
Veremos ahora representaciones conjuntas de individuos y variables, en las que se muestran las contribuciones de los individuos a las componentes principales y el peso de cada variable en las componentes principales.
Variables y observaciones en la primera y segunda componente principal
fviz_pca(PCA,
alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE,
legend.title="Distancia")+theme_bw()
Variables y observaciones en la primera y tercera componente principal
fviz_pca(PCA,axes=c(1,3),
alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE,
legend.title="Distancia")+theme_bw()
Variables y observaciones en la segunda y tercera componente principal
fviz_pca(PCA,axes=c(2,3),
alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE,
legend.title="Distancia")+theme_bw()
Mostramos las coordenadas de los datos originales tipificados en el nuevo sistema de referencia.
PCA$x
## PC1 PC2 PC3 PC4 PC5
## áfrica sur -1.1879510 -1.57174870 0.45044485 0.37814927 1.42283283
## argelia -2.3462462 -1.88106620 -0.81081919 0.20145035 0.45280276
## argentina 1.1979236 -1.53833675 0.20077839 0.89727380 -0.38850796
## australia 3.2068509 -0.89415852 0.83940147 1.06527240 0.55028348
## brasil -0.8909634 -1.34735175 0.89893299 0.55508246 -0.13891942
## canada 2.7566836 -1.13584507 1.63380445 0.83558450 -0.40111280
## chile 0.3748013 -2.00974285 -0.30168062 0.13027061 -1.15642746
## china -2.1745978 2.61356042 0.96542176 0.79447745 -1.51186129
## coreasur -0.2395061 0.22577834 -0.02200379 -0.29949146 -0.46326227
## egipto -2.8332447 -0.89472527 -1.74869384 -0.02886629 0.06798311
## españa 1.9771834 -0.23696615 -1.33663589 0.04438475 -0.47649830
## filipina -2.1452831 0.20351988 1.27403374 -0.86565659 -0.56838811
## francia 2.6414111 -0.01330041 -0.62271201 0.03170101 -0.08010665
## hungría 1.8369158 1.71813286 -0.93757556 0.47001001 0.46001652
## india -3.7665089 1.25869208 1.25543964 -1.00446829 0.51111537
## indonesia -3.5122382 -0.02386264 1.26366487 0.10062250 0.42626703
## iran -2.1106915 -0.95775804 -0.92365774 0.22851570 -0.20529497
## israel 2.5185220 -1.13346621 -0.22444498 -0.43374508 -0.19339430
## italia 1.8793472 0.33682331 -0.08587639 -0.53435523 -0.44626390
## japón 2.4042017 0.77519469 1.67063103 -2.16353425 -0.37124266
## libano 0.3908113 -0.93862865 -1.80755715 -2.54094256 -0.39063809
## marruecos -2.6111083 -0.63494215 -0.84370677 -0.04560012 0.32314717
## méxico -0.4292694 -1.88248069 1.06058617 0.18898306 -0.18870059
## nigeria -4.0691890 0.51312721 1.11817945 0.17003627 0.74676184
## pakistan -3.7752500 -0.29164848 -0.45582169 -0.43747966 0.81514836
## polonia 1.4051850 1.79382646 -0.34557393 0.76060587 0.06958040
## rd alemania 3.0006719 1.80429843 -0.54847315 0.41154507 0.79574301
## reino unido 3.3366283 0.19074891 0.77266251 -1.61255080 0.44377012
## rf alemania 3.6113690 0.93736385 -0.47800299 -1.52413606 1.06149382
## rumania 0.7068599 1.72531630 -0.57683789 0.83431076 0.06419030
## turquia -2.2589795 1.18954704 -2.39282243 0.64594436 -0.65460455
## urss 2.0656447 0.85194470 0.06708247 2.12196985 0.75589667
## usa 2.1968825 -0.80816808 0.55718509 0.86676472 -0.76556968
## vietnam -3.1568661 2.05632216 0.43464714 -0.24212838 -0.56623978
## PC6 PC7 PC8 PC9 PC10
## áfrica sur -0.27199158 -0.623595597 -0.09144981 0.34753493 0.36296233
## argelia -0.21767407 -0.113018266 0.67613362 0.31890906 -0.22213535
## argentina -0.77144470 0.807228821 -0.15129705 -0.40998101 0.13438946
## australia 0.38518213 -0.811803953 0.25022260 -0.16193690 -0.38067934
## brasil 0.02175312 0.175309634 -0.23014686 -0.55123723 0.07607427
## canada 1.26624990 0.214600937 -0.41474660 0.39540041 -0.09652726
## chile -0.18275397 -0.754674607 -0.13156616 -0.14262580 0.15878588
## china 0.21360710 -0.315617302 0.16287419 0.09008596 -0.04968074
## coreasur 1.11284552 -0.035491221 0.39806886 -0.23092496 0.26947888
## egipto 0.02250445 -0.310866567 -0.08354177 -0.02901969 -0.19432866
## españa 0.82678614 0.777742199 0.48167714 -0.10099465 0.01181882
## filipina -0.15028934 0.004493093 0.49597736 0.08712967 0.12559592
## francia 0.13854922 -0.260637363 0.10989111 0.31541253 -0.16761311
## hungría 0.88238111 0.833538559 0.01574888 0.25974597 0.17058796
## india -0.20212478 0.579046678 -0.03383508 -0.25659043 -0.21257778
## indonesia 0.13318440 0.115415733 -0.06815793 0.33036118 -0.08581324
## iran 0.12404861 -0.318863831 0.50032561 -0.17821369 0.21184100
## israel 0.10225569 0.802635989 -0.30878270 -0.14148837 -0.06317868
## italia -1.29071596 0.829905483 0.37519670 0.17010251 -0.10240096
## japón -0.68185453 -0.375862213 -0.00357082 0.01804949 -0.07676834
## libano -0.69319566 0.064675345 -0.38259552 0.27888287 0.13462981
## marruecos 0.86093321 0.077375097 0.12651960 -0.09999321 -0.03204079
## méxico -0.56822183 -0.038134015 0.30275495 -0.16956075 -0.16756015
## nigeria 0.22834537 0.182754443 -0.44626412 -0.22970715 0.15560703
## pakistan -0.13809459 0.282765202 -0.24412230 0.23118350 -0.06426054
## polonia -0.49281049 -0.275363943 0.12995559 -0.09238346 0.01584053
## rd alemania -0.84486979 -0.720954096 -0.03279271 -0.22607825 0.14754470
## reino unido 0.60953258 -0.592216125 -0.31921611 -0.18569892 -0.00228887
## rf alemania 0.72445007 0.104964263 0.11641759 -0.19173155 -0.05768685
## rumania -0.20141636 -0.406726190 0.36289963 0.11225930 0.12367151
## turquia 0.03651040 -0.430409682 -0.82424809 -0.19586549 -0.26359945
## urss -0.97941596 0.614086633 -0.12622830 0.03093806 -0.01039851
## usa -0.07839738 -0.061084782 -0.51133189 0.47419412 0.12224111
## vietnam 0.07615199 -0.021218358 -0.10076961 0.13384194 0.02846942
## PC11
## áfrica sur 0.0009993838
## argelia -0.1452388044
## argentina -0.0562911808
## australia 0.2406669888
## brasil -0.2189263852
## canada -0.1070909000
## chile 0.1209991593
## china 0.0175549674
## coreasur 0.1594077183
## egipto -0.0573909279
## españa 0.1049016670
## filipina 0.2772966071
## francia 0.0703157933
## hungría -0.0509025203
## india -0.1058064933
## indonesia 0.0175770992
## iran -0.2075043694
## israel 0.0268765029
## italia -0.0920630334
## japón -0.2402706122
## libano 0.2033940439
## marruecos -0.2011568416
## méxico 0.1298877117
## nigeria 0.1771866008
## pakistan 0.0721410202
## polonia 0.0187492876
## rd alemania -0.0453764779
## reino unido -0.0604243295
## rf alemania 0.0292461854
## rumania -0.1884743097
## turquia 0.0158359250
## urss 0.1909211317
## usa -0.1407301850
## vietnam 0.0436895774
Como comprobamos anteriormente las correlaciones son significativamente distintas de 0, por lo que tiene sentido el Análisis Factorial.
Podemos visualizar las correlaciones:
##
## Attaching package: 'polycor'
## The following object is masked from 'package:psych':
##
## polyserial
data_fa <- data_norm
poly_cor <- hetcor(data_fa)$correlations
ggcorrplot(poly_cor, type = "lower", hc.order = T)
Ahora, compararemos las salidas de modelos de 3 factores con el método del factor principal y con el de máxima verosimilitud.
modelo1<-fa(poly_cor,
nfactors = 3,
rotate = "none",
fm="mle") # modelo máxima verosimilitud
modelo2<-fa(poly_cor,
nfactors = 3,
rotate = "none",
fm="pa") # modelo factor principal
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
Ante la advertencia decidimos probar el método de mínimo residuo
modelo2<-fa(poly_cor,
nfactors = 3,
rotate = "none",
fm="minres") # modelo mínimo residuo
Finalmente probamos otro método de mínimo residuo
modelo2<-fa(poly_cor,
nfactors = 3,
rotate = "none",
fm="old.min") # modelo mínimo residuo
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
Comparamos las comunalidades:
sort(modelo1$communality,decreasing = T)->c1
sort(modelo2$communality,decreasing = T)->c2
head(cbind(c1,c2))
## c1 c2
## ZPAGRICU 0.9896670 0.9862799
## ZTMINFAN 0.9768999 0.9560850
## ZESPVIDA 0.9632148 0.9192056
## ZTENERGI 0.9577942 0.8752384
## ZPSERVI 0.9286098 0.8726959
## ZPOBURB 0.9120259 0.7414501
Comparamos las unicidades, es decir, la proporción de varianza que no ha sido explicada por el factor (1-comunalidad):
sort(modelo1$uniquenesses,decreasing = T)->u1
sort(modelo2$uniquenesses,decreasing = T)->u2
head(cbind(u1,u2))
## u1 u2
## ZTEJERCI 0.97397969 0.9449885
## ZPOBDENS 0.90544077 0.4806882
## ZTLIBROP 0.41684561 0.4225841
## ZTMEDICO 0.31612567 0.4110643
## ZTPOBACT 0.23641177 0.2962771
## ZPOBURB 0.08797406 0.2585499
Determinaremos el número óptimo de factores por método Scree plot (Cattel 1966) y Análisis Paralelo (Horn 1965):
scree(poly_cor)
fa.parallel(poly_cor, n.obs = nrow(data_fa), fa = "fa", fm = "mle")
## Parallel analysis suggests that the number of factors = 2 and the number of components = NA
El número óptimo de factores es 2.
Estimamos el modelo factorial con 2 factores implementando una rotación tipo varimax para buscar una interpretación más simple.
modelo_varimax<-fa(poly_cor,nfactors = 2,rotate = "varimax", fa="mle")
Mostramos la matriz de pesos factorial rotada
print(modelo_varimax$loadings,cut=0)
##
## Loadings:
## MR1 MR2
## ZPOBDENS 0.023 0.285
## ZTMINFAN -0.756 -0.574
## ZESPVIDA 0.769 0.509
## ZPOBURB 0.964 0.057
## ZTMEDICO 0.631 0.524
## ZPAGRICU -0.985 -0.051
## ZPSERVI 0.941 -0.245
## ZTLIBROP 0.647 0.425
## ZTEJERCI 0.000 0.115
## ZTPOBACT 0.135 0.839
## ZTENERGI 0.603 0.575
##
## MR1 MR2
## SS loadings 5.143 2.239
## Proportion Var 0.468 0.204
## Cumulative Var 0.468 0.671
Veamos una representación más interpretable
fa.diagram(modelo_varimax)
Vemos que el primer factor está asociado con ZPAGRICU, ZPOBURB,ZPSERVI, ZESPVIDA, ZTMINFAN, ZTLIBROP, ZTMEDICO y ZTENERGI, mientras que el segundo factor está asociado con ZTPOBACT.
Veamos con el test de hipótesis que contrasta si el número de factores es suficiente.
factanal(data_fa,factors=2, rotation="none")
##
## Call:
## factanal(x = data_fa, factors = 2, rotation = "none")
##
## Uniquenesses:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## 0.908 0.020 0.053 0.071 0.402 0.051 0.085 0.446
## ZTEJERCI ZTPOBACT ZTENERGI
## 0.987 0.341 0.453
##
## Loadings:
## Factor1 Factor2
## ZPOBDENS 0.147 0.266
## ZTMINFAN -0.960 -0.244
## ZESPVIDA 0.953 0.195
## ZPOBURB 0.889 -0.373
## ZTMEDICO 0.763 0.128
## ZPAGRICU -0.887 0.403
## ZPSERVI 0.758 -0.584
## ZTLIBROP 0.743
## ZTEJERCI 0.111
## ZTPOBACT 0.460 0.669
## ZTENERGI 0.732 0.104
##
## Factor1 Factor2
## SS loadings 5.884 1.298
## Proportion Var 0.535 0.118
## Cumulative Var 0.535 0.653
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 63.09 on 34 degrees of freedom.
## The p-value is 0.00176
Vemos p valor muy bajo, se rechaza la hipótesis nula, 2 no son suficientes.
Estimamos el modelo factorial con 3 factores
Veamos una representación más interpretable
modelo_varimax<-fa(poly_cor,nfactors = 3,rotate = "varimax", fa="mle")
fa.diagram(modelo_varimax)
Vemos que el primer factor está asociado con ZPAGRICU, ZPOBURB,ZPSERVI, ZESPVIDA, ZTMINFAN, ZTLIBROP y ZTMEDICO, mientras que el segundo factor está asociado con ZTPOBACT y ZTENERGI, el tercero con ZTEJERCI.
Veamos con el test de hipótesis que contrasta si el número de factores es suficiente.
factanal(data_fa,factors=3, rotation="none")
##
## Call:
## factanal(x = data_fa, factors = 3, rotation = "none")
##
## Uniquenesses:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## 0.905 0.023 0.037 0.088 0.316 0.010 0.071 0.417
## ZTEJERCI ZTPOBACT ZTENERGI
## 0.974 0.236 0.042
##
## Loadings:
## Factor1 Factor2 Factor3
## ZPOBDENS 0.291
## ZTMINFAN -0.893 -0.395 0.154
## ZESPVIDA 0.888 0.354 -0.221
## ZPOBURB 0.930 -0.199
## ZTMEDICO 0.760 0.269 0.185
## ZPAGRICU -0.961 0.255
## ZPSERVI 0.831 -0.455 -0.178
## ZTLIBROP 0.738 0.180
## ZTEJERCI 0.119 0.105
## ZTPOBACT 0.357 0.775 0.186
## ZTENERGI 0.770 0.281 0.535
##
## Factor1 Factor2 Factor3
## SS loadings 5.915 1.477 0.487
## Proportion Var 0.538 0.134 0.044
## Cumulative Var 0.538 0.672 0.716
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 26.21 on 25 degrees of freedom.
## The p-value is 0.397
No se rechaza la hipótesis nula de que sean suficientes.
Utilizaremos funciones del paquete MVN para contrastar la normalidad multivariante. Veamos si tenemos outliers:
library(MVN)
outliers <- mvn(data = datos_pca , mvnTest = "hz", multivariateOutlierMethod = "quan")
Se detectan 11 outliers en las observaciones. Veamos tests de Royston y Henze-Zirkler:
royston_test <- mvn(data = datos_pca , mvnTest = "royston", multivariatePlot = "qq")
royston_test$multivariateNormality
hz_test <- mvn(data = datos_pca , mvnTest = "hz")
hz_test$multivariateNormality
El test de Royston nos indica rechazar la hip. nula, no tendríamos normalidad multivariante. El test de Henze-Zirkler no nos sugiere rechzar la hip. nula de normalidad multivariante, aunque da p-valor bajo.
# install.packages("tidyverse")
# install.packages("cluster")
# install.packages("factoextra")
# Cargamos los paquetes indicados
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble 3.1.4 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1
## ✓ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x ggplot2::%+%() masks psych::%+%()
## x ggplot2::alpha() masks psych::alpha()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(cluster)
library(factoextra)
data_clust <- datos_pca
Matriz de distancias:
distance<- get_dist(data_clust)
fviz_dist(distance, gradient = list(low ="#00AFBB", mid = "white", high = "#FC4E07"))
Aplicamos clustering con K-Medias:
k2 <- kmeans(data_clust, centers = 2, nstart = 25)
k3 <- kmeans(data_clust, centers = 3, nstart = 25)
k4 <- kmeans(data_clust, centers = 4, nstart = 25)
k5 <- kmeans(data_clust, centers = 5, nstart = 25)
k6 <- kmeans(data_clust, centers = 6, nstart = 25)
# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = data_clust) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = data_clust) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = data_clust) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = data_clust) + ggtitle("k = 5")
p5 <- fviz_cluster(k6, geom = "point", data = data_clust) + ggtitle("k = 6")
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
grid.arrange(p1, p2, p3, p4, p5, nrow = 2)
Buscando el número óptimo de clusters:
set.seed(123)
fviz_nbclust(data_clust, kmeans, method = "wss")
Método de Silhouette
fviz_nbclust(data_clust, kmeans, method = "silhouette")
Método estadístico de brecha (GAP)
gap_stat <- clusGap(data_clust, FUN = kmeans, nstart = 25, K.max = 10, B = 50)
fviz_gap_stat(gap_stat)
Dos de los métodos sugieren que K=2 es el número de clusters más adecuado.
set.seed(123)
final <- kmeans(data_clust, 2, nstart = 25)
print(final)
## K-means clustering with 2 clusters of sizes 18, 16
##
## Cluster means:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## 1 0.2043031 -0.7853941 0.7769316 0.7240592 0.7671309 -0.7448802 0.5956415
## 2 -0.2298410 0.8835684 -0.8740481 -0.8145666 -0.8630222 0.8379902 -0.6700966
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## 1 0.6366902 0.1499519 0.4474332 0.6725565
## 2 -0.7162765 -0.1686959 -0.5033623 -0.7566261
##
## Clustering vector:
## áfrica sur argelia argentina australia brasil canada
## 2 2 1 1 2 1
## chile china coreasur egipto españa filipina
## 1 2 2 2 1 2
## francia hungría india indonesia iran israel
## 1 1 2 2 2 1
## italia japón libano marruecos méxico nigeria
## 1 1 1 2 2 2
## pakistan polonia rd alemania reino unido rf alemania rumania
## 2 1 1 1 1 1
## turquia urss usa vietnam
## 2 1 1 2
##
## Within cluster sum of squares by cluster:
## [1] 107.32705 88.71239
## (between_SS / total_SS = 46.0 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"
# Visualizamos los resultados
fviz_cluster(final, data = data_clust)
Finalmente, mostramos medias de las variables a nivel de cluster.
as.data.frame(data_clust) %>%
mutate(Cluster = final$cluster) %>%
group_by(Cluster) %>%
summarise_all("mean")